function [LE1, LE2] = LEs2D(f, u, k, a, b, x0, q0, Ntotal, Ndiscard)
% LEs2D 计算二维离散系统的两个李雅普诺夫指数
% 输入:
%   f         - 迭代函数句柄，形式为 [x_next, q_next] = f(x, q, u, k, a, b)
%               f 可以为 nhmm、mlm 等
%   u, k, a, b- 系统参数
%   x0, q0   - 初始状态
%   Ntotal   - 总迭代步数
%   Ndiscard - 舍弃的前期步数(不计入 LE 累计)
% 输出:
%   LE1, LE2 - 系统的两个李雅普诺夫指数（按从大到小排序）
%
% 说明:
%   本函数采用有限差分法近似求雅可比矩阵，然后利用 QR 正交化
%   逐步累积对数范数，最后计算李雅普诺夫指数。

    if Ndiscard >= Ntotal
        error('Ndiscard 必须小于 Ntotal！');
    end

    % 初始状态
    x = x0;
    q = q0;
    
    % 初始化正交向量（2D）
    v1 = [1; 0];
    v2 = [0; 1];
    
    sumLog1 = 0;
    sumLog2 = 0;
    
    for n = 1:Ntotal
        % 1. 计算雅可比矩阵：采用有限差分法
        J = jacobian2D(f, u, k, a, b, x, q);
        
        % 2. 对正交向量做线性变换
        v1 = J * v1;
        v2 = J * v2;
        
        % 3. 正交化（QR方法手动实现）
        normv1 = norm(v1);
        v1 = v1 / (normv1 + 1e-30);
        
        proj = dot(v2, v1);
        v2 = v2 - proj * v1;
        normv2 = norm(v2);
        v2 = v2 / (normv2 + 1e-30);
        
        % 4. 更新系统状态
        [x, q] = f(x, q, u, k, a, b);
        
        % 5. 舍弃初始Ndiscard步后累积对数
        if n > Ndiscard
            sumLog1 = sumLog1 + log(normv1 + 1e-30);
            sumLog2 = sumLog2 + log(normv2 + 1e-30);
        end
    end
    
    Ncompute = Ntotal - Ndiscard;
    lambda1 = sumLog1 / Ncompute;
    lambda2 = sumLog2 / Ncompute;
    
    % 排序：通常 lambda1 >= lambda2
    LE = sort([lambda1, lambda2], 'descend');
    LE1 = LE(1);
    LE2 = LE(2);
end

function J = jacobian2D(f, u, k, a, b, x, q)
% 利用有限差分法近似计算二维系统在 (x,q) 处的雅可比矩阵
%   f 为迭代函数句柄：[x_next, q_next] = f(x, q, u, k, a, b)
%   h 为有限差分步长
%
% 输出:
%   J 为 2×2 矩阵，第一列为对 x 的偏导，第二列为对 q 的偏导

    h = 1e-8;  % 小步长
    % 当前 f 值
    [fx, fq] = f(x, q, u, k, a, b);
    
    % 对 x 方向求偏导
    [fx_x, fq_x] = f(x + h, q, u, k, a, b);
    dfdx = ([fx_x; fq_x] - [fx; fq]) / h;
    
    % 对 q 方向求偏导
    [fx_q, fq_q] = f(x, q + h, u, k, a, b);
    dfdq = ([fx_q; fq_q] - [fx; fq]) / h;
    
    J = [dfdx, dfdq];  % 拼接为2×2矩阵
end